Projet Chamois
1 Chargement des librairies
library(tidyverse)
library(corrplot)
library(lmerTest)
library(ade4)
library(splines)
library(car)
library(plotly)
library(DT)
library(Hmisc)
library(kableExtra)
library(knitr)2 Import et description du jeu de données
2.1 Import des données
2.2 Description des données
2.2.1 Résumé des données
## 'data.frame': 1328 obs. of 7 variables:
## $ id : Factor w/ 217 levels "101","105","106",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ year : int 1998 1999 2000 2001 2002 2003 2004 2005 2006 2007 ...
## $ fec : int 1 1 1 1 1 1 1 0 0 0 ...
## $ coh : int 1995 1995 1995 1995 1995 1995 1995 1995 1995 1995 ...
## $ anmark: int 1998 1998 1998 1998 1998 1998 1998 1998 1998 1998 ...
## $ pds : num NA NA NA NA NA NA NA NA NA NA ...
## $ ydth : int 2008 2008 2008 2008 2008 2008 2008 2008 2008 2008 ...
## cham
##
## 7 Variables 1328 Observations
## --------------------------------------------------------------------------------
## id
## n missing distinct
## 1328 0 217
##
## lowest : 101 105 106 107 108, highest: 82 87 9 93 R1
## --------------------------------------------------------------------------------
## year
## n missing distinct Info Mean Gmd .05 .10
## 1328 0 27 0.998 2006 6.831 1995 1997
## .25 .50 .75 .90 .95
## 2001 2006 2010 2014 2015
##
## lowest : 1991 1992 1993 1994 1995, highest: 2013 2014 2015 2016 2017
## --------------------------------------------------------------------------------
## fec
## n missing distinct Info Sum Mean Gmd
## 1328 0 2 0.716 806 0.6069 0.4775
##
## --------------------------------------------------------------------------------
## coh
## n missing distinct Info Mean Gmd .05 .10
## 1328 0 33 0.997 1996 7.75 1985 1987
## .25 .50 .75 .90 .95
## 1991 1997 2001 2005 2007
##
## lowest : 1977 1978 1980 1982 1983, highest: 2007 2009 2010 2011 2014
## --------------------------------------------------------------------------------
## anmark
## n missing distinct Info Mean Gmd .05 .10
## 1328 0 24 0.996 2002 6.288 1993 1994
## .25 .50 .75 .90 .95
## 1998 2002 2006 2009 2011
##
## lowest : 1991 1992 1993 1994 1995, highest: 2010 2011 2012 2014 2015
## --------------------------------------------------------------------------------
## pds
## n missing distinct Info Mean Gmd .05 .10
## 1100 228 92 0.999 19.89 5.25 11.5 12.0
## .25 .50 .75 .90 .95
## 16.9 21.1 23.3 25.0 26.0
##
## lowest : 7.8 10.5 11.0 11.1 11.3, highest: 26.5 26.8 27.0 28.3 28.4
## --------------------------------------------------------------------------------
## ydth
## n missing distinct Info Mean Gmd .05 .10
## 920 408 22 0.977 2006 4.908 1998 2000
## .25 .50 .75 .90 .95
## 2003 2007 2008 2012 2014
##
## lowest : 1994 1996 1997 1998 1999, highest: 2012 2013 2014 2015 2016
## --------------------------------------------------------------------------------
2.2.2 Histogramme années d’observation
2.2.3 Histogramme individus suivis
2.2.4 Présentation des données
Le data frame est constitué de 7 variables et 1328 observations. Chaque observation correspond à l’information de fécondité associée à une femelle chamois et relative à une année donnée. Le jeu de données résume les suivis réalisés entre 1991 et 2017 sur 27 années. D’après l’histogramme des années d’observation, les années entre 2004 et 2006 sont les années pour lesquelles le nombre de chamois suivis a été le plus important (cf histogramme années d’observation). 217 femelles chamois ont été suivies au total. Le nombre d’années de suivi varie selon les femelles entre 1 et 16 années (cf histogramme individus suivis).
2.3 Prétraitement du jeu de données
2.3.1 Elimination des données aberrantes
Les chamois observés après leur mort ou avant leur naissance sont retirés du jeu de donnees.
cham <- cham %>%
filter(year<=ydth | is.na(cham$ydth)) %>%
filter(year>=coh)2.3.2 Création des variables âge (age) et longévité (long)
cham2 <- cham %>%
summarise(cham, age= year-coh, long=ydth-coh, agemark=anmark-coh)3 Question 1 : Lien fécondité annuelle et âge des femelles
3.1 Représentation graphique des données
3.1.1 Représentation par classe d’âge
3.1.2 Représentation sans grouper par classe d’âge
3.1.2.1 Utilisation de la fonction jitter
3.1.2.2 Utilisation de la fonction geom_count
Graphiquement, une augmentation de l’âge des chamois semble engendrer une diminution de la fécondité des chamois avec une inflexion de la courbe observée autour de 10 ans.
3.2 Analyse statistique du lien entre fécondité annuelle et l’âge des femelles
3.2.1 Modèles de régression lineaire generalisé avec effets aléatoires
3.2.1.1 Premier modèle testé
Le premier modèle appliqué est un modèle glm qui utilise la fonction de lien binomial afin de prendre en compte le fait que la variable réponse soit une variable binomiale. La variable “id” est désignée comme variable aléatoire pour prendre en compte le fait que les observations sont répetées sur les mêmes individus sur plusieurs années.
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: fec ~ age + (1 | id)
## Data: cham2
##
## AIC BIC logLik deviance df.resid
## 1761.9 1777.4 -877.9 1755.9 1314
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.6385 -1.1378 0.6862 0.7921 1.0431
##
## Random effects:
## Groups Name Variance Std.Dev.
## id (Intercept) 0.2088 0.4569
## Number of obs: 1317, groups: id, 217
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.58495 0.15766 3.710 0.000207 ***
## age -0.01789 0.01565 -1.143 0.253095
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## age -0.903
| Chisq | Df | Pr(>Chisq) | |
|---|---|---|---|
| age | 1.30614 | 1 | 0.2530947 |
Interpretation des coefficients:
## [1] 1.336377
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.58494620 0.15766366 3.710089 0.0002071865
## age -0.01788695 0.01565097 -1.142865 0.2530947451
## [1] 0.9822721
L’AIC de ce modèle est de 1736. Avec ce modele, la dispersion calculée comme le ratio variance/df est de 1.3 donc il n’y a pas de surdispersion observee. On calcule en utilisant la fonction inverse de logit le coefficient qui permet d’exprimer la fécondité annuelle (en prenant en compte la fonction de lien) en fonction de l’âge. Il est 1.8047875% moins vraisemblable que les chamois aient un petit lorsque leur âge augmente d’un an mais la p value est de 0.25 donc l’impact de l’âge via ce premier modèle n’est pas significatif.
3.2.1.2 Second modèle testé
Un modèle quadratique est testé par la suite pour prendre en compte la tendance de la ligne de régression observée sur les graphiques qui présente une inflexion autour de 10 ans.
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00435823 (tol = 0.002, component 1)
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model is nearly unidentifiable: very large eigenvalue
## - Rescale variables?
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: fec ~ age + I(age^2) + (1 | id)
## Data: cham2
##
## AIC BIC logLik deviance df.resid
## 1674.2 1694.9 -833.1 1666.2 1313
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.1557 -1.0113 0.5808 0.7269 2.8480
##
## Random effects:
## Groups Name Variance Std.Dev.
## id (Intercept) 0.2407 0.4906
## Number of obs: 1317, groups: id, 217
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.076697 0.338088 -6.142 8.12e-10 ***
## age 0.642350 0.077437 8.295 < 2e-16 ***
## I(age^2) -0.033970 0.003999 -8.495 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) age
## age -0.948
## I(age^2) 0.875 -0.977
## optimizer (Nelder_Mead) convergence code: 0 (OK)
## Model failed to converge with max|grad| = 0.00435823 (tol = 0.002, component 1)
## Model is nearly unidentifiable: very large eigenvalue
## - Rescale variables?
La variable age est centrée normée car le modèle n’arrive pas à converger.
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: fec ~ age_scale + I(age_scale^2) + (1 | id)
## Data: cham2
##
## AIC BIC logLik deviance df.resid
## 1674.2 1694.9 -833.1 1666.2 1313
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.1557 -1.0113 0.5808 0.7269 2.8481
##
## Random effects:
## Groups Name Variance Std.Dev.
## id (Intercept) 0.2407 0.4906
## Number of obs: 1317, groups: id, 217
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.95576 0.09371 10.199 <2e-16 ***
## age_scale 0.09431 0.06654 1.417 0.156
## I(age_scale^2) -0.53218 0.06267 -8.492 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) ag_scl
## age_scale 0.146
## I(ag_scl^2) -0.654 -0.168
Interpretation des coefficients:
## [1] 1.26885
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.95576463 0.09371212 10.198944 2.004403e-24
## age_scale 0.09430895 0.06653992 1.417329 1.563868e-01
## I(age_scale^2) -0.53217960 0.06267034 -8.491730 2.035830e-17
L’AIC de ce modèle est de 1674. Avec ce modèle, la dispersion calculée est de 1.3 donc il n’y a pas de surdispersion observée. L’AIC de ce modèle 2 quadratique < l’AIC du modèle 1 linéaire donc le modèle quadratique est plus adapté comme attendu graphiquement. Une observation des coefficients associés aux termes âge et âge^2 indique que le terme “âge” n’est pas significatif dans la prédiction de la variable réponse alors que la p value associée au terme “âge^2” < 0.01. La fonction carré est donc testée.
3.2.1.3 Troisième modèle testé
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: fec ~ I(age_scale^2) + (1 | id)
## Data: cham2
##
## AIC BIC logLik deviance df.resid
## 1674.2 1689.7 -834.1 1668.2 1314
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.1199 -1.0427 0.5836 0.7341 3.0118
##
## Random effects:
## Groups Name Variance Std.Dev.
## id (Intercept) 0.2536 0.5036
## Number of obs: 1317, groups: id, 217
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.93801 0.09350 10.03 <2e-16 ***
## I(age_scale^2) -0.51897 0.06283 -8.26 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## I(ag_scl^2) -0.647
Interpretation des coefficients:
## [1] 1.269406
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.9380093 0.09350033 10.032150 1.100933e-23
## I(age_scale^2) -0.5189717 0.06282985 -8.259955 1.457278e-16
## [1] 0.5951322
L’AIC de ce modèle est de 1674. Avec ce modèle, la dispersion calculée est 1.3 donc il n’y a pas de surdispersion observée. Il est 1.8047875% moins vraisemblable que les chamois aient un petit lorsque leur âge au carré augmente d’une unité (p value < 0.05).
3.2.1.4 4ème modèle testé
On rajoute la variable “year” comme variable aléatoire pour prendre également en compte le fait que les mêmes individus ont été suivis sur les mêmes années.
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: fec ~ I(age_scale^2) + (1 | id) + (1 | year)
## Data: cham2
##
## AIC BIC logLik deviance df.resid
## 1640.9 1661.6 -816.4 1632.9 1313
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.2022 -0.8931 0.5182 0.6924 5.6258
##
## Random effects:
## Groups Name Variance Std.Dev.
## id (Intercept) 0.2863 0.5351
## year (Intercept) 0.3040 0.5514
## Number of obs: 1317, groups: id, 217; year, 27
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.01678 0.14911 6.819 9.17e-12 ***
## I(age_scale^2) -0.56448 0.06598 -8.555 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## I(ag_scl^2) -0.431
Interpretation des coefficients:
## [1] 1.243717
L’AIC de ce modèle est de 1641. Avec ce modèle, la dispersion calculée est de 1.2 donc il n’y a pas de surdispersion observée.
3.2.2 Conclusions
anova(glm1, glm1q, glm1c, glm1d)| npar | AIC | BIC | logLik | deviance | Chisq | Df | Pr(>Chisq) | |
|---|---|---|---|---|---|---|---|---|
| glm1 | 3 | 1761.863 | 1777.412 | -877.9314 | 1755.863 | NA | NA | NA |
| glm1c | 3 | 1674.160 | 1689.710 | -834.0802 | 1668.160 | 87.702339 | 0 | NA |
| glm1q | 4 | 1674.163 | 1694.895 | -833.0813 | 1666.163 | 1.997873 | 1 | 0.1575201 |
| glm1d | 4 | 1640.892 | 1661.624 | -816.4458 | 1632.892 | 33.270918 | 0 | NA |
Le modèle glm1d présente le plus faible AIC. La variable “âge” présente un effet significatif sur la fécondité annuelle via ce modèle ce qui n’est pas surprenant car graphiquement la ligne de régression présentait une courbe avec une diminution de la fécondité pour des âges élevés.
4 Variation de la fécondité annuelle en fonction du temps
4.1 Représentation graphique des données
4.1.1 Représentation graphique par année
4.1.2 Représentation graphique sans grouper par annee
Graphiquement, la variable âge ne semble pas avoir d’effet sur la
fécondité annuelle des chamois malgré le fait que l’âge moyen de la
population augmente sensiblement avec les années.
4.2 Analyse statistique du lien entre fécondité annuelle et années
4.2.1 Modèles de régression lineaire generalisé avec effets aléatoires
4.2.1.1 Premier modèle
Le premier modèle appliqué est un modèle glm qui utilise la fonction de lien binomial afin de prendre en compte le fait que la variable réponse soit une variable binomiale. La variable “id” est désignée comme variable aléatoire pour prendre en compte le fait que les observations sont répetées sur les mêmes individus sur plusieurs années.
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.28972 (tol = 0.002, component 1)
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model is nearly unidentifiable: very large eigenvalue
## - Rescale variables?;Model is nearly unidentifiable: large eigenvalue ratio
## - Rescale variables?
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: fec ~ year_scale + (1 | id)
## Data: cham2
##
## AIC BIC logLik deviance df.resid
## 1763.1 1778.7 -878.6 1757.1 1314
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.6594 -1.1429 0.6860 0.7915 1.0701
##
## Random effects:
## Groups Name Variance Std.Dev.
## id (Intercept) 0.2138 0.4623
## Number of obs: 1317, groups: id, 217
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.42246 0.06787 6.225 4.83e-10 ***
## year_scale -0.01480 0.06540 -0.226 0.821
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## year_scale 0.002
Interpretation des coefficients:
## [1] 1.337139
L’AIC de ce modèle est de 1763. Avec ce modèle, la dispersion calculée comme le ratio variance/df est de 1.3 donc il n’y a pas de surdispersion observée. D’après la p-value > 0.82, il n’y a pas d’effets significatifs de la variable “year” sur la fécondité annuelle comme supposé préalablement par les représentations graphiques. La variable “âge” n’ayant pas d’effets significatifs sur la fécondité annuelle des chamois, il n’y a pas d’utilité à appliquer des modèles avec effets additifs ou multiplicatifs en utilisant les variables âge et année.
5 Lien entre fecondite totale et longevite des animaux
5.1 Représentation graphique des données
5.1.1 Représentation sans prendre en compte le nombre d’annees de suivi
cham_id = cham2 %>%
group_by(id) %>%
dplyr::summarise(feconditetotale= sum(fec), long=long, pds=pds, coh=coh, anneetot=(ydth-min(year)+1), MinAn=min(year), MaxAn=max(year), AgeMark=min(age)) %>%
unique()%>%
drop_na(long)plot9 <-ggplot(cham_id, aes(x=long, y=feconditetotale)) +
geom_count() +
labs(title = "Fecondite totale en fonction de la longevite",x="Longevite", y="Fecondite totale") +
theme(plot.title = element_text(hjust = 0.5)) +
geom_smooth(); plot9
plot10 <- ggplot(cham_id, aes(x=long, y=feconditetotale)) +
geom_jitter(width = 0.25, height = 0.25)+
labs(title = "Fecondite totale en fonction de la longevite",x="Longevite", y="Fecondite totale") +
theme(plot.title = element_text(hjust = 0.5)) +
geom_smooth(); plot10
### Prise en compte du biais apporté par le nombre d’années de suivi
Tous les chamois n’ont pas été marqués au même age et donc n’ont pas été suivis le même nombre d’annee comme le montre l’histogramme ci-dessous.
plot11 <- ggplot(cham_id, aes(x=anneetot)) +
geom_histogram(bins = 30)+
labs(title = "Nombre annees de suivi des chamois",x="Nombre d'annees") +
theme(plot.title = element_text(hjust = 0.5)); plot11plot(jitter(cham_id$long), jitter(cham_id$anneetot))
Le nombre d’annees de suivi n’est donc pas égal à la longévité des
individus comme illustré par le graphique ci-dessous.
plot5 <- ggplot(data = cham_id, aes(x = anneetot,y=long))+
labs(title = "Correlation entre le nombre d'annees de suivi et la longevite",x="Annees de suivi", y="Longevite") +
theme(plot.title = element_text(hjust = 0.5)) +
geom_point()+
geom_smooth(); plot5
On vérifie dans un premier temps, l’impact du nombre d’annees de suivi
sur la fecondite totale des chamois
plot10bis <- ggplot(cham_id, aes(x=anneetot, y=feconditetotale)) +
geom_point(stat="identity") +
labs(title = "Fecondite totale en fonction du nombre d'annee de suivi",x="Nombre annees de suivi", y="Fecondite totale") +
theme(plot.title = element_text(hjust = 0.5)) +
geom_smooth(); plot10bis
Comme attendu, la fecondite totale augmente avec le nombre d’annees de
suivi. Or il est difficile de savoir si cette augmentation est du au
fait que l’individu a été marqué précocement ou si l’individu a vécu
plus longtemps.
Pour pouvoir répondre à la question initiale, qui consiste à vérifier s’il y a un lien entre la fécondité annuelle et la longévité, il faut pouvoir comparer des individus suivis sur le même nombre d’année et si possible sur le maximum d’années possibles.
plot11 <- ggplot(cham_id, aes(x=AgeMark)) +
geom_histogram(bins=50)+
labs(title = "Age de marquage des individus",x="Age") +
theme(plot.title = element_text(hjust = 0.5)); plot11
En s’appuyant sur le graphique ci-dessous, l’analyse statistique va être
réalisé sur un échantillon de la population suivi à savoir tous les
individus marqués à 3 ans qui constitue un échantillon > 30 individus
et peut permettre d’étudier l’impact de la longévité sur la fécondité
totale.
On sait que la variable “année” n’a pas d’impact sur la fécondité annuelle des chamois donc le faut que les chamois aient été suivis pendant des périodes différentes n’a pas d’impact. Il faut par contre vérifier au préalable que le fait de sélectionner la fécondité de la population des individus marqués à 3 ans est représentative de la population totale de chamois.
r <- tapply(cham2$fec, cham2$agemark, mean)
r## 0 1 2 3 4 5 6 7
## 0.5784314 0.6170213 0.6818182 0.6699029 0.6346154 0.5725806 0.6608696 0.6477273
## 8 9 10 11 12 13 14 15
## 0.6071429 0.5740741 0.5573770 0.6065574 0.5483871 0.5454545 0.5000000 1.0000000
## 16 17 21
## 0.2000000 0.0000000 0.0000000
Age3 <- cham2[cham2$agemark==3,]
Ageautre <- cham2[cham2$agemark!=3,]
mean(Age3$fec)## [1] 0.6699029
mean(Ageautre$fec)## [1] 0.6004942
#Voir avec Karim quel test5.1.2 Création du sous-échantillon pour répondre à la question
cham3 <- cham_id %>%
filter(AgeMark==3)5.1.3 Représentation graphique
plot9 <-ggplot(cham3, aes(x=long, y=feconditetotale)) +
geom_point() +
labs(title = "Fecondite totale en fonction de la longevite",x="Longevite", y="Fecondite totale") +
theme(plot.title = element_text(hjust = 0.5)) +
geom_smooth(); plot95.2 Analyse statistique du lien entre la fécondité annuelle et la longevite
5.2.1 Tests de modèles de régression lineaire generalise avec effets aléatoires
5.2.1.1 First
Plusieurs modèles ont été testés => voir avec Karim
lm1 <-lm(feconditetotale~long, data = cham3)
summary(lm1)##
## Call:
## lm(formula = feconditetotale ~ long, data = cham3)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.9408 -1.1635 0.2819 1.0436 3.5826
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2.52025 0.64088 -3.932 0.000275 ***
## long 0.74610 0.06406 11.648 1.86e-15 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.759 on 47 degrees of freedom
## Multiple R-squared: 0.7427, Adjusted R-squared: 0.7372
## F-statistic: 135.7 on 1 and 47 DF, p-value: 1.863e-15
plot(lm1)glmp <- glm(data=cham3, feconditetotale~long, family="poisson")
summary(glmp)##
## Call:
## glm(formula = feconditetotale ~ long, family = "poisson", data = cham3)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.3726 -0.5124 -0.1843 0.6590 1.3798
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.33487 0.22434 -1.493 0.136
## long 0.17122 0.01803 9.494 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 133.056 on 48 degrees of freedom
## Residual deviance: 37.754 on 47 degrees of freedom
## AIC: 186.32
##
## Number of Fisher Scoring iterations: 4
lm2 <- lm(feconditetotale~log(long), data = cham3)
summary(lm2)##
## Call:
## lm(formula = feconditetotale ~ log(long), data = cham3)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.3986 -1.5727 0.0741 1.2620 5.0012
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -7.3397 1.3263 -5.534 1.36e-06 ***
## log(long) 5.5322 0.6116 9.045 7.38e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.095 on 47 degrees of freedom
## Multiple R-squared: 0.6351, Adjusted R-squared: 0.6274
## F-statistic: 81.81 on 1 and 47 DF, p-value: 7.377e-12
plot(lm2)6 Lien entre fecondite annuelle et longevite des animaux
6.1 Représentation graphique des données
plot12 <- ggplot(cham2, aes(x=long, y=fec)) +
geom_count() +
labs(title = "Fécondité annuelle en fonction de la longevite",x="Longevite", y="Fécondité annuelle")+
geom_smooth()+
theme(plot.title = element_text(hjust = 0.5)); plot12plot13 <- ggplot(cham2, aes(x=long, y=fec)) +
geom_jitter(width = 0.55, height=0) +
labs(title = "Fécondité annuelle en fonction de la longevite",x="Longevite", y="Fécondité annuelle")+
geom_smooth()+
theme(plot.title = element_text(hjust = 0.5)); plot13L’allure concave des lignes de régression illustre une augmentation de la fécondité annuelle avec la longévité jusqu’à atteindre un maximum autour de 13-14 ans puis une diminution de la fécondité annuelle.
6.2 Analyse statistique du lien entre fecondite annuelle et longevite des femelles
6.2.1 Modèles de régression lineaire generalise avec effets aléatoires
6.2.1.1 First
Le premier modele applique est un modele glm qui utilise la fonction de lien binomiale afin de prendre en compte le fait que la variable réponse soit une variable binomiale. La variable “id” est désignée comme variable aléatoire pour prendre en compte le fait que les observations sont repetees sur les mêmes individus sur plusieurs annees.
glm11 <- glmer(fec ~ long + (1| id),data=cham2, family = binomial)
long_scale <- scale(x = cham2$long,center = TRUE,scale = TRUE)
glm11 <- glmer(fec ~ long_scale + (1| id),data=cham2, family = binomial)
summary(glm11)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: fec ~ long_scale + (1 | id)
## Data: cham2
##
## AIC BIC logLik deviance df.resid
## 1218.0 1232.5 -606.0 1212.0 906
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.7916 -1.0634 0.6517 0.7850 1.2047
##
## Random effects:
## Groups Name Variance Std.Dev.
## id (Intercept) 0.3199 0.5656
## Number of obs: 909, groups: id, 163
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.38681 0.08681 4.456 8.36e-06 ***
## long_scale 0.04500 0.08512 0.529 0.597
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## long_scale 0.095
Anova(glm11)| Chisq | Df | Pr(>Chisq) | |
|---|---|---|---|
| long_scale | 0.2793972 | 1 | 0.5970966 |
Interpretation des coefficients:
dispersion_glm1 <- 1212/906;dispersion_glm1 #Surdispersion = deviance/df## [1] 1.337748
L’AIC de ce modèle = 1218. Avec ce modele, la dispersion calculée comme le ratio variance/df = 1.3 donc il n’y a pas de surdispersion observee. Avec ce modele, la p value associé à l’impact de la varibale “longevite” sur la fecondite annuelle est de 0.6 donc l’effet de la longevite sur la variable reponse n’est pas significatif.
6.2.1.2 Second
On applique un modèle quadratique pour prendre en compte la tendance la ligne de régression observee sur les graphiques qui présente une inflexion autour de 13-14 ans.
glm11 <- glmer(fec ~ long_scale + I(long_scale^2)+ (1| id),data=cham2, family = binomial)
summary(glm11)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: fec ~ long_scale + I(long_scale^2) + (1 | id)
## Data: cham2
##
## AIC BIC logLik deviance df.resid
## 1209.3 1228.6 -600.7 1201.3 905
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.7429 -1.0969 0.6542 0.7821 1.4323
##
## Random effects:
## Groups Name Variance Std.Dev.
## id (Intercept) 0.2492 0.4992
## Number of obs: 909, groups: id, 163
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.60021 0.10574 5.676 1.38e-08 ***
## long_scale 0.03059 0.08331 0.367 0.71352
## I(long_scale^2) -0.20676 0.06313 -3.275 0.00106 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) lng_sc
## long_scale 0.026
## I(lng_sc^2) -0.612 0.065
Anova(glm11)| Chisq | Df | Pr(>Chisq) | |
|---|---|---|---|
| long_scale | 0.1347913 | 1 | 0.7135150 |
| I(long_scale^2) | 10.7277777 | 1 | 0.0010554 |
Interpretation des coefficients:
dispersion_glm1 <- 1201/905;dispersion_glm1 #Surdispersion = deviance/df## [1] 1.327072
L’AIC de ce modèle glm1q = 1209. Avec le modele glm1q, la dispersion calculée = 1.3 donc il n’y a pas de surdispersion observee. L’AIC de ce modèle glm1q < l’AIC du modèle glm1 donc le modèle glm1q est plus adapté. Une observation des coefficients indique que le terme “x” n’est pas significatif dans la prediction de la variable reponse alors que la p value associée au terme “x2” < 0.01. La fonction carré est donc testée.
6.2.1.3 Third
glm11 <- glmer(fec ~ I(long_scale^2)+ (1| id),data=cham2, family = binomial)
summary(glm11)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: fec ~ I(long_scale^2) + (1 | id)
## Data: cham2
##
## AIC BIC logLik deviance df.resid
## 1207.5 1221.9 -600.7 1201.5 906
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.7388 -1.0973 0.6529 0.7774 1.3841
##
## Random effects:
## Groups Name Variance Std.Dev.
## id (Intercept) 0.2565 0.5064
## Number of obs: 909, groups: id, 163
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.59931 0.10618 5.644 1.66e-08 ***
## I(long_scale^2) -0.20842 0.06328 -3.294 0.000989 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## I(lng_sc^2) -0.616
Anova(glm11)| Chisq | Df | Pr(>Chisq) | |
|---|---|---|---|
| I(long_scale^2) | 10.84755 | 1 | 0.0009893 |
Interpretation des coefficients:
dispersion_glm1c <- 1201/906;dispersion_glm1c #Surdispersion = deviance/df## [1] 1.325607
summary(glm11)$coefficients## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.5993126 0.10618061 5.644275 1.658785e-08
## I(long_scale^2) -0.2084191 0.06328076 -3.293563 9.892629e-04
exp(summary(glm11)$coefficients[[2]])## [1] 0.8118667
Coeffb <- ((1/exp(summary(glm11)$coefficients[2]))-1)*100L’AIC de ce modèle glm1c = 1674. Avec le modele glm1c, la dispersion calculée = 1.3 donc il n’y a pas de surdispersion observee. Il est 23.1729325% moins vraisemblable que les chamois aient un petit lorsque leur longevite au carré augmente d’une unité (p value < 0.05).
Segmentation => voir avec Karim
cham_long1 <- cham2[cham2$long<15,]
cham_long2 <- cham2[cham2$long>=15,]
plot16 <- ggplot(cham_long1, aes(x=long, y=fec)) +
geom_count() +
labs(title = "Fécondité en fonction de la longevite",x="Longevite", y="Fécondité")+
geom_smooth()+
theme(plot.title = element_text(hjust = 0.5)); plot16## Warning: Removed 408 rows containing non-finite values (`stat_sum()`).
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'
## Warning: Removed 408 rows containing non-finite values (`stat_smooth()`).
glm12 <- glmer(fec ~ long + (1| id),data=cham_long1, family = binomial)
summary(glm12)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: fec ~ long + (1 | id)
## Data: cham_long1
##
## AIC BIC logLik deviance df.resid
## 789.0 802.2 -391.5 783.0 585
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.6864 -1.0755 0.6737 0.8230 1.1503
##
## Random effects:
## Groups Name Variance Std.Dev.
## id (Intercept) 0.1568 0.3959
## Number of obs: 588, groups: id, 120
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.66433 0.39918 -1.664 0.09606 .
## long 0.09835 0.03665 2.684 0.00728 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## long -0.971
Anova(glm12)| Chisq | Df | Pr(>Chisq) | |
|---|---|---|---|
| long | 7.201696 | 1 | 0.0072835 |
plot17 <- ggplot(cham_long2, aes(x=long, y=fec)) +
geom_count() +
labs(title = "Fécondité en fonction de la longevite",x="Longevite", y="Fécondité")+
geom_smooth()+
theme(plot.title = element_text(hjust = 0.5)); plot17## Warning: Removed 408 rows containing non-finite values (`stat_sum()`).
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'
## Warning: Removed 408 rows containing non-finite values (`stat_smooth()`).
glm13 <- glmer(fec ~ long + (1| id),data=cham_long2, family = binomial)
summary(glm13)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: fec ~ long + (1 | id)
## Data: cham_long2
##
## AIC BIC logLik deviance df.resid
## 418.9 430.3 -206.5 412.9 318
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.9648 -1.0836 0.5871 0.7880 1.2014
##
## Random effects:
## Groups Name Variance Std.Dev.
## id (Intercept) 0.3949 0.6284
## Number of obs: 321, groups: id, 43
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 4.35885 1.43312 3.041 0.00235 **
## long -0.23441 0.08423 -2.783 0.00539 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## long -0.994
Anova(glm13)| Chisq | Df | Pr(>Chisq) | |
|---|---|---|---|
| long | 7.744764 | 1 | 0.0053868 |
7 Lien entre fecondite totale et poids
7.1 Représentation graphique des données
7.1.1 Vérification de la comparabilite des poids selon les ages de capture
cham_pds <- cham_id %>%
drop_na(pds)plot18 <- ggplot(cham_pds, aes(x=AgeMark, y=pds)) +
geom_point(color="blue", fill=rgb(0.1,0.4,0.5,0.7)) +
labs(title = "Poids selon l'age de capture",x="Age de capture", y="Poids mesuré")+
theme(plot.title = element_text(hjust = 0.5))+
geom_smooth();plot18#Fiare une boucle for pour exclure outlier
(boxplot(cham_pds$pds))## $stats
## [,1]
## [1,] 16.0
## [2,] 20.0
## [3,] 22.0
## [4,] 23.5
## [5,] 27.0
##
## $n
## [1] 131
##
## $conf
## [,1]
## [1,] 21.51684
## [2,] 22.48316
##
## $out
## [1] 14.3 12.8 7.8 14.5 12.5 14.0 14.2 11.9 12.5 14.0 11.1 11.3 11.1 14.5 11.5
## [16] 10.5
##
## $group
## [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
##
## $names
## [1] "1"
out <- boxplot.stats(cham_pds$pds)$out
max(out)## [1] 14.5
cham_pds <- cham_pds%>%
filter(pds>max(out))
hist(cham_pds$AgeMark)plot(cham_pds$pds, cham_pds$feconditetotale)plot(cham_pds$AgeMark, cham_pds$feconditetotale)
Reflechir à la section des donnes avant analyse stat pour fecondite
totale
7.2 Analyse statistique du lien entre fecondite annuelle/longevite et poids des femelles
7.2.1 Modèles de régression lineaire generalise avec effets aléatoires
7.2.1.1 First
lm14 <- lm(feconditetotale ~ pds,data=cham_pds)
summary(lm14)##
## Call:
## lm(formula = feconditetotale ~ pds, data = cham_pds)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.205 -2.086 -0.979 1.790 9.910
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.28615 2.26237 1.011 0.314
## pds 0.03828 0.10097 0.379 0.705
##
## Residual standard error: 2.526 on 113 degrees of freedom
## Multiple R-squared: 0.00127, Adjusted R-squared: -0.007568
## F-statistic: 0.1437 on 1 and 113 DF, p-value: 0.7053
lm15 <- lm(long ~ pds,data=cham_pds)
summary(lm15)##
## Call:
## lm(formula = long ~ pds, data = cham_pds)
##
## Residuals:
## Min 1Q Median 3Q Max
## -8.1325 -2.3945 0.0646 2.3182 7.8675
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.4877 3.0358 -0.490 0.625
## pds 0.6010 0.1355 4.435 2.14e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.389 on 113 degrees of freedom
## Multiple R-squared: 0.1483, Adjusted R-squared: 0.1407
## F-statistic: 19.67 on 1 and 113 DF, p-value: 2.143e-05
7.2.2 Représentation graphique de la fecondite totale en fonction du poids
plot19 <- ggplot(cham_pds, aes(x=pds, y=feconditetotale)) +
geom_point(color="blue", fill=rgb(0.1,0.4,0.5,0.7)) +
geom_smooth();plot197.2.3 Représentation graphique de la longevite en fonction du poids
plot20 <- ggplot(cham_pds, aes(x=pds, y=long)) +
geom_point(color="blue", fill=rgb(0.1,0.4,0.5,0.7)) +
geom_smooth();plot20
Longevite significatif sur le poids mais pas de la fecondite totale.